{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {
     "autoexec": {
      "startup": false,
      "wait_interval": 0
     }
    },
    "colab_type": "code",
    "id": "mmKBfDbTa_4M"
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.ndimage "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {
     "autoexec": {
      "startup": false,
      "wait_interval": 0
     }
    },
    "colab_type": "code",
    "id": "-r22-U0ZnrRz"
   },
   "outputs": [],
   "source": [
    "# neighbour_code_to_normals is a lookup table.\n",
    "# For every binary neighbour code \n",
    "# (2x2x2 neighbourhood = 8 neighbours = 8 bits = 256 codes) \n",
    "# it contains the surface normals of the triangles (called \"surfel\" for \n",
    "# \"surface element\" in the following). The length of the normal \n",
    "# vector encodes the surfel area.\n",
    "#\n",
    "# created by compute_surface_area_lookup_table.ipynb using the \n",
    "# marching_cube algorithm, see e.g. https://en.wikipedia.org/wiki/Marching_cubes\n",
    "#\n",
    "neighbour_code_to_normals = [\n",
    "  [[0,0,0]],\n",
    "  [[0.125,0.125,0.125]],\n",
    "  [[-0.125,-0.125,0.125]],\n",
    "  [[-0.25,-0.25,0.0],[0.25,0.25,-0.0]],\n",
    "  [[0.125,-0.125,0.125]],\n",
    "  [[-0.25,-0.0,-0.25],[0.25,0.0,0.25]],\n",
    "  [[0.125,-0.125,0.125],[-0.125,-0.125,0.125]],\n",
    "  [[0.5,0.0,-0.0],[0.25,0.25,0.25],[0.125,0.125,0.125]],\n",
    "  [[-0.125,0.125,0.125]],\n",
    "  [[0.125,0.125,0.125],[-0.125,0.125,0.125]],\n",
    "  [[-0.25,0.0,0.25],[-0.25,0.0,0.25]],\n",
    "  [[0.5,0.0,0.0],[-0.25,-0.25,0.25],[-0.125,-0.125,0.125]],\n",
    "  [[0.25,-0.25,0.0],[0.25,-0.25,0.0]],\n",
    "  [[0.5,0.0,0.0],[0.25,-0.25,0.25],[-0.125,0.125,-0.125]],\n",
    "  [[-0.5,0.0,0.0],[-0.25,0.25,0.25],[-0.125,0.125,0.125]],\n",
    "  [[0.5,0.0,0.0],[0.5,0.0,0.0]],\n",
    "  [[0.125,-0.125,-0.125]],\n",
    "  [[0.0,-0.25,-0.25],[0.0,0.25,0.25]],\n",
    "  [[-0.125,-0.125,0.125],[0.125,-0.125,-0.125]],\n",
    "  [[0.0,-0.5,0.0],[0.25,0.25,0.25],[0.125,0.125,0.125]],\n",
    "  [[0.125,-0.125,0.125],[0.125,-0.125,-0.125]],\n",
    "  [[0.0,0.0,-0.5],[0.25,0.25,0.25],[-0.125,-0.125,-0.125]],\n",
    "  [[-0.125,-0.125,0.125],[0.125,-0.125,0.125],[0.125,-0.125,-0.125]],\n",
    "  [[-0.125,-0.125,-0.125],[-0.25,-0.25,-0.25],[0.25,0.25,0.25],[0.125,0.125,0.125]],\n",
    "  [[-0.125,0.125,0.125],[0.125,-0.125,-0.125]],\n",
    "  [[0.0,-0.25,-0.25],[0.0,0.25,0.25],[-0.125,0.125,0.125]],\n",
    "  [[-0.25,0.0,0.25],[-0.25,0.0,0.25],[0.125,-0.125,-0.125]],\n",
    "  [[0.125,0.125,0.125],[0.375,0.375,0.375],[0.0,-0.25,0.25],[-0.25,0.0,0.25]],\n",
    "  [[0.125,-0.125,-0.125],[0.25,-0.25,0.0],[0.25,-0.25,0.0]],\n",
    "  [[0.375,0.375,0.375],[0.0,0.25,-0.25],[-0.125,-0.125,-0.125],[-0.25,0.25,0.0]],\n",
    "  [[-0.5,0.0,0.0],[-0.125,-0.125,-0.125],[-0.25,-0.25,-0.25],[0.125,0.125,0.125]],\n",
    "  [[-0.5,0.0,0.0],[-0.125,-0.125,-0.125],[-0.25,-0.25,-0.25]],\n",
    "  [[0.125,-0.125,0.125]],\n",
    "  [[0.125,0.125,0.125],[0.125,-0.125,0.125]],\n",
    "  [[0.0,-0.25,0.25],[0.0,0.25,-0.25]],\n",
    "  [[0.0,-0.5,0.0],[0.125,0.125,-0.125],[0.25,0.25,-0.25]],\n",
    "  [[0.125,-0.125,0.125],[0.125,-0.125,0.125]],\n",
    "  [[0.125,-0.125,0.125],[-0.25,-0.0,-0.25],[0.25,0.0,0.25]],\n",
    "  [[0.0,-0.25,0.25],[0.0,0.25,-0.25],[0.125,-0.125,0.125]],\n",
    "  [[-0.375,-0.375,0.375],[-0.0,0.25,0.25],[0.125,0.125,-0.125],[-0.25,-0.0,-0.25]],\n",
    "  [[-0.125,0.125,0.125],[0.125,-0.125,0.125]],\n",
    "  [[0.125,0.125,0.125],[0.125,-0.125,0.125],[-0.125,0.125,0.125]],\n",
    "  [[-0.0,0.0,0.5],[-0.25,-0.25,0.25],[-0.125,-0.125,0.125]],\n",
    "  [[0.25,0.25,-0.25],[0.25,0.25,-0.25],[0.125,0.125,-0.125],[-0.125,-0.125,0.125]],\n",
    "  [[0.125,-0.125,0.125],[0.25,-0.25,0.0],[0.25,-0.25,0.0]],\n",
    "  [[0.5,0.0,0.0],[0.25,-0.25,0.25],[-0.125,0.125,-0.125],[0.125,-0.125,0.125]],\n",
    "  [[0.0,0.25,-0.25],[0.375,-0.375,-0.375],[-0.125,0.125,0.125],[0.25,0.25,0.0]],\n",
    "  [[-0.5,0.0,0.0],[-0.25,-0.25,0.25],[-0.125,-0.125,0.125]],\n",
    "  [[0.25,-0.25,0.0],[-0.25,0.25,0.0]],\n",
    "  [[0.0,0.5,0.0],[-0.25,0.25,0.25],[0.125,-0.125,-0.125]],\n",
    "  [[0.0,0.5,0.0],[0.125,-0.125,0.125],[-0.25,0.25,-0.25]],\n",
    "  [[0.0,0.5,0.0],[0.0,-0.5,0.0]],\n",
    "  [[0.25,-0.25,0.0],[-0.25,0.25,0.0],[0.125,-0.125,0.125]],\n",
    "  [[-0.375,-0.375,-0.375],[-0.25,0.0,0.25],[-0.125,-0.125,-0.125],[-0.25,0.25,0.0]],\n",
    "  [[0.125,0.125,0.125],[0.0,-0.5,0.0],[-0.25,-0.25,-0.25],[-0.125,-0.125,-0.125]],\n",
    "  [[0.0,-0.5,0.0],[-0.25,-0.25,-0.25],[-0.125,-0.125,-0.125]],\n",
    "  [[-0.125,0.125,0.125],[0.25,-0.25,0.0],[-0.25,0.25,0.0]],\n",
    "  [[0.0,0.5,0.0],[0.25,0.25,-0.25],[-0.125,-0.125,0.125],[-0.125,-0.125,0.125]],\n",
    "  [[-0.375,0.375,-0.375],[-0.25,-0.25,0.0],[-0.125,0.125,-0.125],[-0.25,0.0,0.25]],\n",
    "  [[0.0,0.5,0.0],[0.25,0.25,-0.25],[-0.125,-0.125,0.125]],\n",
    "  [[0.25,-0.25,0.0],[-0.25,0.25,0.0],[0.25,-0.25,0.0],[0.25,-0.25,0.0]],\n",
    "  [[-0.25,-0.25,0.0],[-0.25,-0.25,0.0],[-0.125,-0.125,0.125]],\n",
    "  [[0.125,0.125,0.125],[-0.25,-0.25,0.0],[-0.25,-0.25,0.0]],\n",
    "  [[-0.25,-0.25,0.0],[-0.25,-0.25,0.0]],\n",
    "  [[-0.125,-0.125,0.125]],\n",
    "  [[0.125,0.125,0.125],[-0.125,-0.125,0.125]],\n",
    "  [[-0.125,-0.125,0.125],[-0.125,-0.125,0.125]],\n",
    "  [[-0.125,-0.125,0.125],[-0.25,-0.25,0.0],[0.25,0.25,-0.0]],\n",
    "  [[0.0,-0.25,0.25],[0.0,-0.25,0.25]],\n",
    "  [[0.0,0.0,0.5],[0.25,-0.25,0.25],[0.125,-0.125,0.125]],\n",
    "  [[0.0,-0.25,0.25],[0.0,-0.25,0.25],[-0.125,-0.125,0.125]],\n",
    "  [[0.375,-0.375,0.375],[0.0,-0.25,-0.25],[-0.125,0.125,-0.125],[0.25,0.25,0.0]],\n",
    "  [[-0.125,-0.125,0.125],[-0.125,0.125,0.125]],\n",
    "  [[0.125,0.125,0.125],[-0.125,-0.125,0.125],[-0.125,0.125,0.125]],\n",
    "  [[-0.125,-0.125,0.125],[-0.25,0.0,0.25],[-0.25,0.0,0.25]],\n",
    "  [[0.5,0.0,0.0],[-0.25,-0.25,0.25],[-0.125,-0.125,0.125],[-0.125,-0.125,0.125]],\n",
    "  [[-0.0,0.5,0.0],[-0.25,0.25,-0.25],[0.125,-0.125,0.125]],\n",
    "  [[-0.25,0.25,-0.25],[-0.25,0.25,-0.25],[-0.125,0.125,-0.125],[-0.125,0.125,-0.125]],\n",
    "  [[-0.25,0.0,-0.25],[0.375,-0.375,-0.375],[0.0,0.25,-0.25],[-0.125,0.125,0.125]],\n",
    "  [[0.5,0.0,0.0],[-0.25,0.25,-0.25],[0.125,-0.125,0.125]],\n",
    "  [[-0.25,0.0,0.25],[0.25,0.0,-0.25]],\n",
    "  [[-0.0,0.0,0.5],[-0.25,0.25,0.25],[-0.125,0.125,0.125]],\n",
    "  [[-0.125,-0.125,0.125],[-0.25,0.0,0.25],[0.25,0.0,-0.25]],\n",
    "  [[-0.25,-0.0,-0.25],[-0.375,0.375,0.375],[-0.25,-0.25,0.0],[-0.125,0.125,0.125]],\n",
    "  [[0.0,0.0,-0.5],[0.25,0.25,-0.25],[-0.125,-0.125,0.125]],\n",
    "  [[-0.0,0.0,0.5],[0.0,0.0,0.5]],\n",
    "  [[0.125,0.125,0.125],[0.125,0.125,0.125],[0.25,0.25,0.25],[0.0,0.0,0.5]],\n",
    "  [[0.125,0.125,0.125],[0.25,0.25,0.25],[0.0,0.0,0.5]],\n",
    "  [[-0.25,0.0,0.25],[0.25,0.0,-0.25],[-0.125,0.125,0.125]],\n",
    "  [[-0.0,0.0,0.5],[0.25,-0.25,0.25],[0.125,-0.125,0.125],[0.125,-0.125,0.125]],\n",
    "  [[-0.25,0.0,0.25],[-0.25,0.0,0.25],[-0.25,0.0,0.25],[0.25,0.0,-0.25]],\n",
    "  [[0.125,-0.125,0.125],[0.25,0.0,0.25],[0.25,0.0,0.25]],\n",
    "  [[0.25,0.0,0.25],[-0.375,-0.375,0.375],[-0.25,0.25,0.0],[-0.125,-0.125,0.125]],\n",
    "  [[-0.0,0.0,0.5],[0.25,-0.25,0.25],[0.125,-0.125,0.125]],\n",
    "  [[0.125,0.125,0.125],[0.25,0.0,0.25],[0.25,0.0,0.25]],\n",
    "  [[0.25,0.0,0.25],[0.25,0.0,0.25]],\n",
    "  [[-0.125,-0.125,0.125],[0.125,-0.125,0.125]],\n",
    "  [[0.125,0.125,0.125],[-0.125,-0.125,0.125],[0.125,-0.125,0.125]],\n",
    "  [[-0.125,-0.125,0.125],[0.0,-0.25,0.25],[0.0,0.25,-0.25]],\n",
    "  [[0.0,-0.5,0.0],[0.125,0.125,-0.125],[0.25,0.25,-0.25],[-0.125,-0.125,0.125]],\n",
    "  [[0.0,-0.25,0.25],[0.0,-0.25,0.25],[0.125,-0.125,0.125]],\n",
    "  [[0.0,0.0,0.5],[0.25,-0.25,0.25],[0.125,-0.125,0.125],[0.125,-0.125,0.125]],\n",
    "  [[0.0,-0.25,0.25],[0.0,-0.25,0.25],[0.0,-0.25,0.25],[0.0,0.25,-0.25]],\n",
    "  [[0.0,0.25,0.25],[0.0,0.25,0.25],[0.125,-0.125,-0.125]],\n",
    "  [[-0.125,0.125,0.125],[0.125,-0.125,0.125],[-0.125,-0.125,0.125]],\n",
    "  [[-0.125,0.125,0.125],[0.125,-0.125,0.125],[-0.125,-0.125,0.125],[0.125,0.125,0.125]],\n",
    "  [[-0.0,0.0,0.5],[-0.25,-0.25,0.25],[-0.125,-0.125,0.125],[-0.125,-0.125,0.125]],\n",
    "  [[0.125,0.125,0.125],[0.125,-0.125,0.125],[0.125,-0.125,-0.125]],\n",
    "  [[-0.0,0.5,0.0],[-0.25,0.25,-0.25],[0.125,-0.125,0.125],[0.125,-0.125,0.125]],\n",
    "  [[0.125,0.125,0.125],[-0.125,-0.125,0.125],[0.125,-0.125,-0.125]],\n",
    "  [[0.0,-0.25,-0.25],[0.0,0.25,0.25],[0.125,0.125,0.125]],\n",
    "  [[0.125,0.125,0.125],[0.125,-0.125,-0.125]],\n",
    "  [[0.5,0.0,-0.0],[0.25,-0.25,-0.25],[0.125,-0.125,-0.125]],\n",
    "  [[-0.25,0.25,0.25],[-0.125,0.125,0.125],[-0.25,0.25,0.25],[0.125,-0.125,-0.125]],\n",
    "  [[0.375,-0.375,0.375],[0.0,0.25,0.25],[-0.125,0.125,-0.125],[-0.25,0.0,0.25]],\n",
    "  [[0.0,-0.5,0.0],[-0.25,0.25,0.25],[-0.125,0.125,0.125]],\n",
    "  [[-0.375,-0.375,0.375],[0.25,-0.25,0.0],[0.0,0.25,0.25],[-0.125,-0.125,0.125]],\n",
    "  [[-0.125,0.125,0.125],[-0.25,0.25,0.25],[0.0,0.0,0.5]],\n",
    "  [[0.125,0.125,0.125],[0.0,0.25,0.25],[0.0,0.25,0.25]],\n",
    "  [[0.0,0.25,0.25],[0.0,0.25,0.25]],\n",
    "  [[0.5,0.0,-0.0],[0.25,0.25,0.25],[0.125,0.125,0.125],[0.125,0.125,0.125]],\n",
    "  [[0.125,-0.125,0.125],[-0.125,-0.125,0.125],[0.125,0.125,0.125]],\n",
    "  [[-0.25,-0.0,-0.25],[0.25,0.0,0.25],[0.125,0.125,0.125]],\n",
    "  [[0.125,0.125,0.125],[0.125,-0.125,0.125]],\n",
    "  [[-0.25,-0.25,0.0],[0.25,0.25,-0.0],[0.125,0.125,0.125]],\n",
    "  [[0.125,0.125,0.125],[-0.125,-0.125,0.125]],\n",
    "  [[0.125,0.125,0.125],[0.125,0.125,0.125]],\n",
    "  [[0.125,0.125,0.125]],\n",
    "  [[0.125,0.125,0.125]],\n",
    "  [[0.125,0.125,0.125],[0.125,0.125,0.125]],\n",
    "  [[0.125,0.125,0.125],[-0.125,-0.125,0.125]],\n",
    "  [[-0.25,-0.25,0.0],[0.25,0.25,-0.0],[0.125,0.125,0.125]],\n",
    "  [[0.125,0.125,0.125],[0.125,-0.125,0.125]],\n",
    "  [[-0.25,-0.0,-0.25],[0.25,0.0,0.25],[0.125,0.125,0.125]],\n",
    "  [[0.125,-0.125,0.125],[-0.125,-0.125,0.125],[0.125,0.125,0.125]],\n",
    "  [[0.5,0.0,-0.0],[0.25,0.25,0.25],[0.125,0.125,0.125],[0.125,0.125,0.125]],\n",
    "  [[0.0,0.25,0.25],[0.0,0.25,0.25]],\n",
    "  [[0.125,0.125,0.125],[0.0,0.25,0.25],[0.0,0.25,0.25]],\n",
    "  [[-0.125,0.125,0.125],[-0.25,0.25,0.25],[0.0,0.0,0.5]],\n",
    "  [[-0.375,-0.375,0.375],[0.25,-0.25,0.0],[0.0,0.25,0.25],[-0.125,-0.125,0.125]],\n",
    "  [[0.0,-0.5,0.0],[-0.25,0.25,0.25],[-0.125,0.125,0.125]],\n",
    "  [[0.375,-0.375,0.375],[0.0,0.25,0.25],[-0.125,0.125,-0.125],[-0.25,0.0,0.25]],\n",
    "  [[-0.25,0.25,0.25],[-0.125,0.125,0.125],[-0.25,0.25,0.25],[0.125,-0.125,-0.125]],\n",
    "  [[0.5,0.0,-0.0],[0.25,-0.25,-0.25],[0.125,-0.125,-0.125]],\n",
    "  [[0.125,0.125,0.125],[0.125,-0.125,-0.125]],\n",
    "  [[0.0,-0.25,-0.25],[0.0,0.25,0.25],[0.125,0.125,0.125]],\n",
    "  [[0.125,0.125,0.125],[-0.125,-0.125,0.125],[0.125,-0.125,-0.125]],\n",
    "  [[-0.0,0.5,0.0],[-0.25,0.25,-0.25],[0.125,-0.125,0.125],[0.125,-0.125,0.125]],\n",
    "  [[0.125,0.125,0.125],[0.125,-0.125,0.125],[0.125,-0.125,-0.125]],\n",
    "  [[-0.0,0.0,0.5],[-0.25,-0.25,0.25],[-0.125,-0.125,0.125],[-0.125,-0.125,0.125]],\n",
    "  [[-0.125,0.125,0.125],[0.125,-0.125,0.125],[-0.125,-0.125,0.125],[0.125,0.125,0.125]],\n",
    "  [[-0.125,0.125,0.125],[0.125,-0.125,0.125],[-0.125,-0.125,0.125]],\n",
    "  [[0.0,0.25,0.25],[0.0,0.25,0.25],[0.125,-0.125,-0.125]],\n",
    "  [[0.0,-0.25,-0.25],[0.0,0.25,0.25],[0.0,0.25,0.25],[0.0,0.25,0.25]],\n",
    "  [[0.0,0.0,0.5],[0.25,-0.25,0.25],[0.125,-0.125,0.125],[0.125,-0.125,0.125]],\n",
    "  [[0.0,-0.25,0.25],[0.0,-0.25,0.25],[0.125,-0.125,0.125]],\n",
    "  [[0.0,-0.5,0.0],[0.125,0.125,-0.125],[0.25,0.25,-0.25],[-0.125,-0.125,0.125]],\n",
    "  [[-0.125,-0.125,0.125],[0.0,-0.25,0.25],[0.0,0.25,-0.25]],\n",
    "  [[0.125,0.125,0.125],[-0.125,-0.125,0.125],[0.125,-0.125,0.125]],\n",
    "  [[-0.125,-0.125,0.125],[0.125,-0.125,0.125]],\n",
    "  [[0.25,0.0,0.25],[0.25,0.0,0.25]],\n",
    "  [[0.125,0.125,0.125],[0.25,0.0,0.25],[0.25,0.0,0.25]],\n",
    "  [[-0.0,0.0,0.5],[0.25,-0.25,0.25],[0.125,-0.125,0.125]],\n",
    "  [[0.25,0.0,0.25],[-0.375,-0.375,0.375],[-0.25,0.25,0.0],[-0.125,-0.125,0.125]],\n",
    "  [[0.125,-0.125,0.125],[0.25,0.0,0.25],[0.25,0.0,0.25]],\n",
    "  [[-0.25,-0.0,-0.25],[0.25,0.0,0.25],[0.25,0.0,0.25],[0.25,0.0,0.25]],\n",
    "  [[-0.0,0.0,0.5],[0.25,-0.25,0.25],[0.125,-0.125,0.125],[0.125,-0.125,0.125]],\n",
    "  [[-0.25,0.0,0.25],[0.25,0.0,-0.25],[-0.125,0.125,0.125]],\n",
    "  [[0.125,0.125,0.125],[0.25,0.25,0.25],[0.0,0.0,0.5]],\n",
    "  [[0.125,0.125,0.125],[0.125,0.125,0.125],[0.25,0.25,0.25],[0.0,0.0,0.5]],\n",
    "  [[-0.0,0.0,0.5],[0.0,0.0,0.5]],\n",
    "  [[0.0,0.0,-0.5],[0.25,0.25,-0.25],[-0.125,-0.125,0.125]],\n",
    "  [[-0.25,-0.0,-0.25],[-0.375,0.375,0.375],[-0.25,-0.25,0.0],[-0.125,0.125,0.125]],\n",
    "  [[-0.125,-0.125,0.125],[-0.25,0.0,0.25],[0.25,0.0,-0.25]],\n",
    "  [[-0.0,0.0,0.5],[-0.25,0.25,0.25],[-0.125,0.125,0.125]],\n",
    "  [[-0.25,0.0,0.25],[0.25,0.0,-0.25]],\n",
    "  [[0.5,0.0,0.0],[-0.25,0.25,-0.25],[0.125,-0.125,0.125]],\n",
    "  [[-0.25,0.0,-0.25],[0.375,-0.375,-0.375],[0.0,0.25,-0.25],[-0.125,0.125,0.125]],\n",
    "  [[-0.25,0.25,-0.25],[-0.25,0.25,-0.25],[-0.125,0.125,-0.125],[-0.125,0.125,-0.125]],\n",
    "  [[-0.0,0.5,0.0],[-0.25,0.25,-0.25],[0.125,-0.125,0.125]],\n",
    "  [[0.5,0.0,0.0],[-0.25,-0.25,0.25],[-0.125,-0.125,0.125],[-0.125,-0.125,0.125]],\n",
    "  [[-0.125,-0.125,0.125],[-0.25,0.0,0.25],[-0.25,0.0,0.25]],\n",
    "  [[0.125,0.125,0.125],[-0.125,-0.125,0.125],[-0.125,0.125,0.125]],\n",
    "  [[-0.125,-0.125,0.125],[-0.125,0.125,0.125]],\n",
    "  [[0.375,-0.375,0.375],[0.0,-0.25,-0.25],[-0.125,0.125,-0.125],[0.25,0.25,0.0]],\n",
    "  [[0.0,-0.25,0.25],[0.0,-0.25,0.25],[-0.125,-0.125,0.125]],\n",
    "  [[0.0,0.0,0.5],[0.25,-0.25,0.25],[0.125,-0.125,0.125]],\n",
    "  [[0.0,-0.25,0.25],[0.0,-0.25,0.25]],\n",
    "  [[-0.125,-0.125,0.125],[-0.25,-0.25,0.0],[0.25,0.25,-0.0]],\n",
    "  [[-0.125,-0.125,0.125],[-0.125,-0.125,0.125]],\n",
    "  [[0.125,0.125,0.125],[-0.125,-0.125,0.125]],\n",
    "  [[-0.125,-0.125,0.125]],\n",
    "  [[-0.25,-0.25,0.0],[-0.25,-0.25,0.0]],\n",
    "  [[0.125,0.125,0.125],[-0.25,-0.25,0.0],[-0.25,-0.25,0.0]],\n",
    "  [[-0.25,-0.25,0.0],[-0.25,-0.25,0.0],[-0.125,-0.125,0.125]],\n",
    "  [[-0.25,-0.25,0.0],[-0.25,-0.25,0.0],[-0.25,-0.25,0.0],[0.25,0.25,-0.0]],\n",
    "  [[0.0,0.5,0.0],[0.25,0.25,-0.25],[-0.125,-0.125,0.125]],\n",
    "  [[-0.375,0.375,-0.375],[-0.25,-0.25,0.0],[-0.125,0.125,-0.125],[-0.25,0.0,0.25]],\n",
    "  [[0.0,0.5,0.0],[0.25,0.25,-0.25],[-0.125,-0.125,0.125],[-0.125,-0.125,0.125]],\n",
    "  [[-0.125,0.125,0.125],[0.25,-0.25,0.0],[-0.25,0.25,0.0]],\n",
    "  [[0.0,-0.5,0.0],[-0.25,-0.25,-0.25],[-0.125,-0.125,-0.125]],\n",
    "  [[0.125,0.125,0.125],[0.0,-0.5,0.0],[-0.25,-0.25,-0.25],[-0.125,-0.125,-0.125]],\n",
    "  [[-0.375,-0.375,-0.375],[-0.25,0.0,0.25],[-0.125,-0.125,-0.125],[-0.25,0.25,0.0]],\n",
    "  [[0.25,-0.25,0.0],[-0.25,0.25,0.0],[0.125,-0.125,0.125]],\n",
    "  [[0.0,0.5,0.0],[0.0,-0.5,0.0]],\n",
    "  [[0.0,0.5,0.0],[0.125,-0.125,0.125],[-0.25,0.25,-0.25]],\n",
    "  [[0.0,0.5,0.0],[-0.25,0.25,0.25],[0.125,-0.125,-0.125]],\n",
    "  [[0.25,-0.25,0.0],[-0.25,0.25,0.0]],\n",
    "  [[-0.5,0.0,0.0],[-0.25,-0.25,0.25],[-0.125,-0.125,0.125]],\n",
    "  [[0.0,0.25,-0.25],[0.375,-0.375,-0.375],[-0.125,0.125,0.125],[0.25,0.25,0.0]],\n",
    "  [[0.5,0.0,0.0],[0.25,-0.25,0.25],[-0.125,0.125,-0.125],[0.125,-0.125,0.125]],\n",
    "  [[0.125,-0.125,0.125],[0.25,-0.25,0.0],[0.25,-0.25,0.0]],\n",
    "  [[0.25,0.25,-0.25],[0.25,0.25,-0.25],[0.125,0.125,-0.125],[-0.125,-0.125,0.125]],\n",
    "  [[-0.0,0.0,0.5],[-0.25,-0.25,0.25],[-0.125,-0.125,0.125]],\n",
    "  [[0.125,0.125,0.125],[0.125,-0.125,0.125],[-0.125,0.125,0.125]],\n",
    "  [[-0.125,0.125,0.125],[0.125,-0.125,0.125]],\n",
    "  [[-0.375,-0.375,0.375],[-0.0,0.25,0.25],[0.125,0.125,-0.125],[-0.25,-0.0,-0.25]],\n",
    "  [[0.0,-0.25,0.25],[0.0,0.25,-0.25],[0.125,-0.125,0.125]],\n",
    "  [[0.125,-0.125,0.125],[-0.25,-0.0,-0.25],[0.25,0.0,0.25]],\n",
    "  [[0.125,-0.125,0.125],[0.125,-0.125,0.125]],\n",
    "  [[0.0,-0.5,0.0],[0.125,0.125,-0.125],[0.25,0.25,-0.25]],\n",
    "  [[0.0,-0.25,0.25],[0.0,0.25,-0.25]],\n",
    "  [[0.125,0.125,0.125],[0.125,-0.125,0.125]],\n",
    "  [[0.125,-0.125,0.125]],\n",
    "  [[-0.5,0.0,0.0],[-0.125,-0.125,-0.125],[-0.25,-0.25,-0.25]],\n",
    "  [[-0.5,0.0,0.0],[-0.125,-0.125,-0.125],[-0.25,-0.25,-0.25],[0.125,0.125,0.125]],\n",
    "  [[0.375,0.375,0.375],[0.0,0.25,-0.25],[-0.125,-0.125,-0.125],[-0.25,0.25,0.0]],\n",
    "  [[0.125,-0.125,-0.125],[0.25,-0.25,0.0],[0.25,-0.25,0.0]],\n",
    "  [[0.125,0.125,0.125],[0.375,0.375,0.375],[0.0,-0.25,0.25],[-0.25,0.0,0.25]],\n",
    "  [[-0.25,0.0,0.25],[-0.25,0.0,0.25],[0.125,-0.125,-0.125]],\n",
    "  [[0.0,-0.25,-0.25],[0.0,0.25,0.25],[-0.125,0.125,0.125]],\n",
    "  [[-0.125,0.125,0.125],[0.125,-0.125,-0.125]],\n",
    "  [[-0.125,-0.125,-0.125],[-0.25,-0.25,-0.25],[0.25,0.25,0.25],[0.125,0.125,0.125]],\n",
    "  [[-0.125,-0.125,0.125],[0.125,-0.125,0.125],[0.125,-0.125,-0.125]],\n",
    "  [[0.0,0.0,-0.5],[0.25,0.25,0.25],[-0.125,-0.125,-0.125]],\n",
    "  [[0.125,-0.125,0.125],[0.125,-0.125,-0.125]],\n",
    "  [[0.0,-0.5,0.0],[0.25,0.25,0.25],[0.125,0.125,0.125]],\n",
    "  [[-0.125,-0.125,0.125],[0.125,-0.125,-0.125]],\n",
    "  [[0.0,-0.25,-0.25],[0.0,0.25,0.25]],\n",
    "  [[0.125,-0.125,-0.125]],\n",
    "  [[0.5,0.0,0.0],[0.5,0.0,0.0]],\n",
    "  [[-0.5,0.0,0.0],[-0.25,0.25,0.25],[-0.125,0.125,0.125]],\n",
    "  [[0.5,0.0,0.0],[0.25,-0.25,0.25],[-0.125,0.125,-0.125]],\n",
    "  [[0.25,-0.25,0.0],[0.25,-0.25,0.0]],\n",
    "  [[0.5,0.0,0.0],[-0.25,-0.25,0.25],[-0.125,-0.125,0.125]],\n",
    "  [[-0.25,0.0,0.25],[-0.25,0.0,0.25]],\n",
    "  [[0.125,0.125,0.125],[-0.125,0.125,0.125]],\n",
    "  [[-0.125,0.125,0.125]],\n",
    "  [[0.5,0.0,-0.0],[0.25,0.25,0.25],[0.125,0.125,0.125]],\n",
    "  [[0.125,-0.125,0.125],[-0.125,-0.125,0.125]],\n",
    "  [[-0.25,-0.0,-0.25],[0.25,0.0,0.25]],\n",
    "  [[0.125,-0.125,0.125]],\n",
    "  [[-0.25,-0.25,0.0],[0.25,0.25,-0.0]],\n",
    "  [[-0.125,-0.125,0.125]],\n",
    "  [[0.125,0.125,0.125]],\n",
    "  [[0,0,0]]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {
     "autoexec": {
      "startup": false,
      "wait_interval": 0
     }
    },
    "colab_type": "code",
    "id": "-kruspp5Y1Ip"
   },
   "outputs": [],
   "source": [
    "def compute_surface_distances(mask_gt, mask_pred, spacing_mm):\n",
    "  \"\"\"Compute closest distances from all surface points to the other surface.\n",
    "\n",
    "  Finds all surface elements \"surfels\" in the ground truth mask `mask_gt` and\n",
    "  the predicted mask `mask_pred`, computes their area in mm^2 and the distance\n",
    "  to the closest point on the other surface. It returns two sorted lists of\n",
    "  distances together with the corresponding surfel areas. If one of the masks\n",
    "  is empty, the corresponding lists are empty and all distances in the other\n",
    "  list are `inf` \n",
    "  \n",
    "  Args:\n",
    "    mask_gt: 3-dim Numpy array of type bool. The ground truth mask.\n",
    "    mask_pred: 3-dim Numpy array of type bool. The predicted mask.\n",
    "    spacing_mm: 3-element list-like structure. Voxel spacing in x0, x1 and x2\n",
    "        direction \n",
    "\n",
    "  Returns:\n",
    "    A dict with \n",
    "    \"distances_gt_to_pred\": 1-dim numpy array of type float. The distances in mm\n",
    "        from all ground truth surface elements to the predicted surface, \n",
    "        sorted from smallest to largest\n",
    "    \"distances_pred_to_gt\": 1-dim numpy array of type float. The distances in mm\n",
    "        from all predicted surface elements to the ground truth surface, \n",
    "        sorted from smallest to largest \n",
    "    \"surfel_areas_gt\": 1-dim numpy array of type float. The area in mm^2 of \n",
    "        the ground truth surface elements in the same order as \n",
    "        distances_gt_to_pred\n",
    "    \"surfel_areas_pred\": 1-dim numpy array of type float. The area in mm^2 of \n",
    "        the predicted surface elements in the same order as \n",
    "        distances_pred_to_gt\n",
    "       \n",
    "  \"\"\"\n",
    "  \n",
    "  # compute the area for all 256 possible surface elements \n",
    "  # (given a 2x2x2 neighbourhood) according to the spacing_mm\n",
    "  neighbour_code_to_surface_area = np.zeros([256])\n",
    "  for code in range(256):\n",
    "    normals = np.array(neighbour_code_to_normals[code])\n",
    "    sum_area = 0\n",
    "    for normal_idx in range(normals.shape[0]):\n",
    "      # normal vector\n",
    "      n = np.zeros([3])\n",
    "      n[0] = normals[normal_idx,0] * spacing_mm[1] * spacing_mm[2]\n",
    "      n[1] = normals[normal_idx,1] * spacing_mm[0] * spacing_mm[2]\n",
    "      n[2] = normals[normal_idx,2] * spacing_mm[0] * spacing_mm[1]\n",
    "      area = np.linalg.norm(n)\n",
    "      sum_area += area\n",
    "    neighbour_code_to_surface_area[code] = sum_area\n",
    "\n",
    "  # compute the bounding box of the masks to trim\n",
    "  # the volume to the smallest possible processing subvolume\n",
    "  mask_all = mask_gt | mask_pred\n",
    "  bbox_min = np.zeros(3, np.int64)\n",
    "  bbox_max = np.zeros(3, np.int64)\n",
    "\n",
    "  # max projection to the x0-axis\n",
    "  proj_0 = np.max(np.max(mask_all, axis=2), axis=1)\n",
    "  idx_nonzero_0 = np.nonzero(proj_0)[0]\n",
    "  if len(idx_nonzero_0) == 0:\n",
    "    return {\"distances_gt_to_pred\":  np.array([]), \n",
    "            \"distances_pred_to_gt\":  np.array([]), \n",
    "            \"surfel_areas_gt\":       np.array([]), \n",
    "            \"surfel_areas_pred\":     np.array([])}\n",
    "    \n",
    "  bbox_min[0] = np.min(idx_nonzero_0)\n",
    "  bbox_max[0] = np.max(idx_nonzero_0)\n",
    "\n",
    "  # max projection to the x1-axis\n",
    "  proj_1 = np.max(np.max(mask_all, axis=2), axis=0)\n",
    "  idx_nonzero_1 = np.nonzero(proj_1)[0]\n",
    "  bbox_min[1] = np.min(idx_nonzero_1)\n",
    "  bbox_max[1] = np.max(idx_nonzero_1)\n",
    "\n",
    "  # max projection to the x2-axis\n",
    "  proj_2 = np.max(np.max(mask_all, axis=1), axis=0)\n",
    "  idx_nonzero_2 = np.nonzero(proj_2)[0]\n",
    "  bbox_min[2] = np.min(idx_nonzero_2)\n",
    "  bbox_max[2] = np.max(idx_nonzero_2)\n",
    "\n",
    "  print(\"bounding box min = {}\".format(bbox_min))\n",
    "  print(\"bounding box max = {}\".format(bbox_max))\n",
    "\n",
    "  # crop the processing subvolume.\n",
    "  # we need to zeropad the cropped region with 1 voxel at the lower, \n",
    "  # the right and the back side. This is required to obtain the \"full\" \n",
    "  # convolution result with the 2x2x2 kernel\n",
    "  cropmask_gt = np.zeros((bbox_max - bbox_min)+2, np.uint8)\n",
    "  cropmask_pred = np.zeros((bbox_max - bbox_min)+2, np.uint8)\n",
    "\n",
    "  cropmask_gt[0:-1, 0:-1, 0:-1] = mask_gt[bbox_min[0]:bbox_max[0]+1,\n",
    "                                          bbox_min[1]:bbox_max[1]+1,\n",
    "                                          bbox_min[2]:bbox_max[2]+1]\n",
    "\n",
    "  cropmask_pred[0:-1, 0:-1, 0:-1] = mask_pred[bbox_min[0]:bbox_max[0]+1,\n",
    "                                              bbox_min[1]:bbox_max[1]+1,\n",
    "                                              bbox_min[2]:bbox_max[2]+1]\n",
    "\n",
    "  # compute the neighbour code (local binary pattern) for each voxel\n",
    "  # the resultsing arrays are spacially shifted by minus half a voxel in each axis.\n",
    "  # i.e. the points are located at the corners of the original voxels\n",
    "  kernel = np.array([[[128,64],\n",
    "                      [32,16]],\n",
    "                     [[8,4],\n",
    "                      [2,1]]])\n",
    "  neighbour_code_map_gt = scipy.ndimage.filters.correlate(cropmask_gt.astype(np.uint8), kernel, mode=\"constant\", cval=0) \n",
    "  neighbour_code_map_pred = scipy.ndimage.filters.correlate(cropmask_pred.astype(np.uint8), kernel, mode=\"constant\", cval=0) \n",
    "\n",
    "  # create masks with the surface voxels\n",
    "  borders_gt   = ((neighbour_code_map_gt != 0) & (neighbour_code_map_gt != 255))\n",
    "  borders_pred = ((neighbour_code_map_pred != 0) & (neighbour_code_map_pred != 255))\n",
    "\n",
    "  # compute the distance transform (closest distance of each voxel to the surface voxels)\n",
    "  if borders_gt.any():\n",
    "    distmap_gt = scipy.ndimage.morphology.distance_transform_edt(~borders_gt, sampling=spacing_mm)\n",
    "  else:\n",
    "    distmap_gt = np.Inf * np.ones(borders_gt.shape)\n",
    "\n",
    "  if borders_pred.any():  \n",
    "    distmap_pred = scipy.ndimage.morphology.distance_transform_edt(~borders_pred, sampling=spacing_mm)\n",
    "  else:\n",
    "    distmap_pred = np.Inf * np.ones(borders_pred.shape)\n",
    "\n",
    "  # compute the area of each surface element\n",
    "  surface_area_map_gt = neighbour_code_to_surface_area[neighbour_code_map_gt]\n",
    "  surface_area_map_pred = neighbour_code_to_surface_area[neighbour_code_map_pred]\n",
    "\n",
    "  # create a list of all surface elements with distance and area\n",
    "  distances_gt_to_pred = distmap_pred[borders_gt]\n",
    "  distances_pred_to_gt = distmap_gt[borders_pred]\n",
    "  surfel_areas_gt   = surface_area_map_gt[borders_gt]\n",
    "  surfel_areas_pred = surface_area_map_pred[borders_pred]\n",
    "\n",
    "  # sort them by distance\n",
    "  if distances_gt_to_pred.shape != (0,):\n",
    "    sorted_surfels_gt = np.array(sorted(zip(distances_gt_to_pred, surfel_areas_gt)))\n",
    "    distances_gt_to_pred = sorted_surfels_gt[:,0]\n",
    "    surfel_areas_gt      = sorted_surfels_gt[:,1]\n",
    "\n",
    "  if distances_pred_to_gt.shape != (0,):\n",
    "    sorted_surfels_pred = np.array(sorted(zip(distances_pred_to_gt, surfel_areas_pred)))\n",
    "    distances_pred_to_gt = sorted_surfels_pred[:,0]\n",
    "    surfel_areas_pred    = sorted_surfels_pred[:,1]\n",
    "\n",
    "\n",
    "  return {\"distances_gt_to_pred\":  distances_gt_to_pred, \n",
    "          \"distances_pred_to_gt\":  distances_pred_to_gt, \n",
    "          \"surfel_areas_gt\":       surfel_areas_gt, \n",
    "          \"surfel_areas_pred\":     surfel_areas_pred}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {
     "autoexec": {
      "startup": false,
      "wait_interval": 0
     }
    },
    "colab_type": "code",
    "id": "y8shBcxlrUtJ"
   },
   "outputs": [],
   "source": [
    "def compute_average_surface_distance(surface_distances):\n",
    "  distances_gt_to_pred = surface_distances[\"distances_gt_to_pred\"]\n",
    "  distances_pred_to_gt = surface_distances[\"distances_pred_to_gt\"]\n",
    "  surfel_areas_gt      = surface_distances[\"surfel_areas_gt\"]\n",
    "  surfel_areas_pred    = surface_distances[\"surfel_areas_pred\"]\n",
    "  average_distance_gt_to_pred = np.sum( distances_gt_to_pred * surfel_areas_gt) / np.sum(surfel_areas_gt)\n",
    "  average_distance_pred_to_gt = np.sum( distances_pred_to_gt * surfel_areas_pred) / np.sum(surfel_areas_pred)\n",
    "  return (average_distance_gt_to_pred, average_distance_pred_to_gt)\n",
    "\n",
    "def compute_robust_hausdorff(surface_distances, percent):\n",
    "  distances_gt_to_pred = surface_distances[\"distances_gt_to_pred\"]\n",
    "  distances_pred_to_gt = surface_distances[\"distances_pred_to_gt\"]\n",
    "  surfel_areas_gt      = surface_distances[\"surfel_areas_gt\"]\n",
    "  surfel_areas_pred    = surface_distances[\"surfel_areas_pred\"]\n",
    "  if len(distances_gt_to_pred) > 0:\n",
    "    surfel_areas_cum_gt   = np.cumsum(surfel_areas_gt) / np.sum(surfel_areas_gt)\n",
    "    idx = np.searchsorted(surfel_areas_cum_gt, percent/100.0)\n",
    "    perc_distance_gt_to_pred = distances_gt_to_pred[min(idx, len(distances_gt_to_pred)-1)]\n",
    "  else:\n",
    "    perc_distance_gt_to_pred = np.Inf\n",
    "    \n",
    "  if len(distances_pred_to_gt) > 0:\n",
    "    surfel_areas_cum_pred = np.cumsum(surfel_areas_pred) / np.sum(surfel_areas_pred)\n",
    "    idx = np.searchsorted(surfel_areas_cum_pred, percent/100.0)\n",
    "    perc_distance_pred_to_gt = distances_pred_to_gt[min(idx, len(distances_pred_to_gt)-1)]\n",
    "  else:\n",
    "    perc_distance_pred_to_gt = np.Inf\n",
    "    \n",
    "  return max( perc_distance_gt_to_pred, perc_distance_pred_to_gt)\n",
    "\n",
    "def compute_surface_overlap_at_tolerance(surface_distances, tolerance_mm):\n",
    "  distances_gt_to_pred = surface_distances[\"distances_gt_to_pred\"]\n",
    "  distances_pred_to_gt = surface_distances[\"distances_pred_to_gt\"]\n",
    "  surfel_areas_gt      = surface_distances[\"surfel_areas_gt\"]\n",
    "  surfel_areas_pred    = surface_distances[\"surfel_areas_pred\"]\n",
    "  rel_overlap_gt   = np.sum(surfel_areas_gt[distances_gt_to_pred <= tolerance_mm]) / np.sum(surfel_areas_gt)\n",
    "  rel_overlap_pred = np.sum(surfel_areas_pred[distances_pred_to_gt <= tolerance_mm]) / np.sum(surfel_areas_pred)\n",
    "  return (rel_overlap_gt, rel_overlap_pred)\n",
    "\n",
    "def compute_surface_dice_at_tolerance(surface_distances, tolerance_mm):\n",
    "  distances_gt_to_pred = surface_distances[\"distances_gt_to_pred\"]\n",
    "  distances_pred_to_gt = surface_distances[\"distances_pred_to_gt\"]\n",
    "  surfel_areas_gt      = surface_distances[\"surfel_areas_gt\"]\n",
    "  surfel_areas_pred    = surface_distances[\"surfel_areas_pred\"]\n",
    "  overlap_gt   = np.sum(surfel_areas_gt[distances_gt_to_pred <= tolerance_mm])\n",
    "  overlap_pred = np.sum(surfel_areas_pred[distances_pred_to_gt <= tolerance_mm])\n",
    "  surface_dice = (overlap_gt + overlap_pred) / (\n",
    "      np.sum(surfel_areas_gt) + np.sum(surfel_areas_pred))\n",
    "  return surface_dice"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {
     "autoexec": {
      "startup": false,
      "wait_interval": 0
     }
    },
    "colab_type": "code",
    "id": "8XgG0aHkHIej"
   },
   "outputs": [],
   "source": [
    "def compute_dice_coefficient(mask_gt, mask_pred):\n",
    "  \"\"\"Compute soerensen-dice coefficient.\n",
    "\n",
    "  compute the soerensen-dice coefficient between the ground truth mask `mask_gt`\n",
    "  and the predicted mask `mask_pred`. \n",
    "  \n",
    "  Args:\n",
    "    mask_gt: 3-dim Numpy array of type bool. The ground truth mask.\n",
    "    mask_pred: 3-dim Numpy array of type bool. The predicted mask.\n",
    "\n",
    "  Returns:\n",
    "    the dice coeffcient as float. If both masks are empty, the result is NaN\n",
    "  \"\"\"\n",
    "  volume_sum = mask_gt.sum() + mask_pred.sum()\n",
    "  if volume_sum == 0:\n",
    "    return np.NaN\n",
    "  volume_intersect = (mask_gt & mask_pred).sum()\n",
    "  return 2*volume_intersect / volume_sum\n",
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "4VLjvdslp1ZW"
   },
   "source": [
    "# Some Simple Tests"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {
     "autoexec": {
      "startup": false,
      "wait_interval": 0
     },
     "base_uri": "https://localhost:8080/",
     "height": 153
    },
    "colab_type": "code",
    "executionInfo": {
     "elapsed": 520,
     "status": "ok",
     "timestamp": 1528999072261,
     "user": {
      "displayName": "Olaf Ronneberger",
      "photoUrl": "//lh5.googleusercontent.com/-Hp2EnPvUOmU/AAAAAAAAAAI/AAAAAAAAAFc/QR20zxrgwZ0/s50-c-k-no/photo.jpg",
      "userId": "110531354013919487725"
     },
     "user_tz": -60
    },
    "id": "Vl25cjRakoSb",
    "outputId": "8b9bec94-7fa4-41cf-d8bf-f1a92db4dda8"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "bounding box min = [50 60 70]\n",
      "bounding box max = [50 60 72]\n",
      "average surface distance: (1.5, 1.5) mm\n",
      "hausdorff (100%):         2.0 mm\n",
      "hausdorff (95%):          2.0 mm\n",
      "surface overlap at 1mm:   (0.5, 0.5)\n",
      "surface dice at 1mm:      0.5\n",
      "volumetric dice:          0.0\n"
     ]
    }
   ],
   "source": [
    "# single pixels, 2mm away\n",
    "mask_gt   = np.zeros((128,128,128), np.uint8)\n",
    "mask_pred = np.zeros((128,128,128), np.uint8)\n",
    "mask_gt[50,60,70] = 1\n",
    "mask_pred[50,60,72] = 1\n",
    "surface_distances = compute_surface_distances(mask_gt, mask_pred, spacing_mm=(3,2,1))\n",
    "print(\"surface dice at 1mm:      {}\".format(compute_surface_dice_at_tolerance(surface_distances, 1)))\n",
    "print(\"volumetric dice:          {}\".format(compute_dice_coefficient(mask_gt, mask_pred)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {
     "autoexec": {
      "startup": false,
      "wait_interval": 0
     },
     "base_uri": "https://localhost:8080/",
     "height": 204
    },
    "colab_type": "code",
    "executionInfo": {
     "elapsed": 763,
     "status": "ok",
     "timestamp": 1528999073141,
     "user": {
      "displayName": "Olaf Ronneberger",
      "photoUrl": "//lh5.googleusercontent.com/-Hp2EnPvUOmU/AAAAAAAAAAI/AAAAAAAAAFc/QR20zxrgwZ0/s50-c-k-no/photo.jpg",
      "userId": "110531354013919487725"
     },
     "user_tz": -60
    },
    "id": "HKp9I52Ik26R",
    "outputId": "f1ecf88c-8114-4139-f5a9-885330249a97"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "bounding box min = [0 0 0]\n",
      "bounding box max = [50 99 99]\n",
      "average surface distance: (0.3224470726766679, 0.3393962930605405) mm\n",
      "hausdorff (100%):         2.0 mm\n",
      "hausdorff (95%):          2.0 mm\n",
      "surface overlap at 1mm:   (0.8420663973547752, 0.8303018534697297)\n",
      "surface dice at 1mm:      0.836145008498\n",
      "volumetric dice:          0.990099009901\n",
      "\n",
      "expected average_distance_gt_to_pred = 1./6 * 2mm = 0.333333333333mm\n",
      "expected volumetric dice: 0.990099009901\n"
     ]
    }
   ],
   "source": [
    "# two cubes. cube 1 is 100x100x100 mm^3 and cube 2 is 102x100x100 mm^3\n",
    "mask_gt   = np.zeros((100,100,100), np.uint8)\n",
    "mask_pred = np.zeros((100,100,100), np.uint8)\n",
    "spacing_mm=(2,1,1)\n",
    "mask_gt[0:50, :, :] = 1\n",
    "mask_pred[0:51, :, :] = 1\n",
    "surface_distances = compute_surface_distances(mask_gt, mask_pred, spacing_mm)\n",
    "print(\"surface dice at 1mm:      {}\".format(compute_surface_dice_at_tolerance(surface_distances, 1)))\n",
    "print(\"volumetric dice:          {}\".format(compute_dice_coefficient(mask_gt, mask_pred)))\n",
    "print(\"\")\n",
    "print(\"expected average_distance_gt_to_pred = 1./6 * 2mm = {}mm\".format(1./6 * 2))  \n",
    "print(\"expected volumetric dice: {}\".format(2.*100*100*100 / (100*100*100 + 102*100*100) ))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {
     "autoexec": {
      "startup": false,
      "wait_interval": 0
     },
     "base_uri": "https://localhost:8080/",
     "height": 204
    },
    "colab_type": "code",
    "executionInfo": {
     "elapsed": 463,
     "status": "ok",
     "timestamp": 1528999073635,
     "user": {
      "displayName": "Olaf Ronneberger",
      "photoUrl": "//lh5.googleusercontent.com/-Hp2EnPvUOmU/AAAAAAAAAAI/AAAAAAAAAFc/QR20zxrgwZ0/s50-c-k-no/photo.jpg",
      "userId": "110531354013919487725"
     },
     "user_tz": -60
    },
    "id": "1EgXUCqNqI3X",
    "outputId": "6f7946b5-1bc3-4cc4-b489-82801e7aad2c"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "bounding box min = [50 60 70]\n",
      "bounding box max = [50 60 70]\n",
      "average surface distance: (inf, nan) mm\n",
      "hausdorff (100%):         inf mm\n",
      "hausdorff (95%):          inf mm\n",
      "surface overlap at 1mm:   (0.0, nan)\n",
      "surface dice at 1mm:      0.0\n",
      "volumetric dice:          0.0\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/lib/python2.7/dist-packages/ipykernel_launcher.py:7: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  import sys\n",
      "/usr/local/lib/python2.7/dist-packages/ipykernel_launcher.py:37: RuntimeWarning: invalid value encountered in double_scalars\n"
     ]
    }
   ],
   "source": [
    "# test empty mask in prediction\n",
    "mask_gt   = np.zeros((128,128,128), np.uint8)\n",
    "mask_pred = np.zeros((128,128,128), np.uint8)\n",
    "mask_gt[50,60,70] = 1\n",
    "#mask_pred[50,60,72] = 1\n",
    "surface_distances = compute_surface_distances(mask_gt, mask_pred, spacing_mm=(3,2,1))\n",
    "print(\"average surface distance: {} mm\".format(compute_average_surface_distance(surface_distances)))\n",
    "print(\"hausdorff (100%):         {} mm\".format(compute_robust_hausdorff(surface_distances, 100)))\n",
    "print(\"hausdorff (95%):          {} mm\".format(compute_robust_hausdorff(surface_distances, 95)))\n",
    "print(\"surface overlap at 1mm:   {}\".format(compute_surface_overlap_at_tolerance(surface_distances, 1)))\n",
    "print(\"surface dice at 1mm:      {}\".format(compute_surface_dice_at_tolerance(surface_distances, 1)))\n",
    "print(\"volumetric dice:          {}\".format(compute_dice_coefficient(mask_gt, mask_pred)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {
     "autoexec": {
      "startup": false,
      "wait_interval": 0
     },
     "base_uri": "https://localhost:8080/",
     "height": 204
    },
    "colab_type": "code",
    "executionInfo": {
     "elapsed": 468,
     "status": "ok",
     "timestamp": 1528999074187,
     "user": {
      "displayName": "Olaf Ronneberger",
      "photoUrl": "//lh5.googleusercontent.com/-Hp2EnPvUOmU/AAAAAAAAAAI/AAAAAAAAAFc/QR20zxrgwZ0/s50-c-k-no/photo.jpg",
      "userId": "110531354013919487725"
     },
     "user_tz": -60
    },
    "id": "TcBi7tcNaK3w",
    "outputId": "a0f488d1-d11a-4a71-f6be-38f8dbaca57b"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "bounding box min = [50 60 72]\n",
      "bounding box max = [50 60 72]\n",
      "average surface distance: (nan, inf) mm\n",
      "hausdorff (100%):         inf mm\n",
      "hausdorff (95%):          inf mm\n",
      "surface overlap at 1mm:   (nan, 0.0)\n",
      "surface dice at 1mm:      0.0\n",
      "volumetric dice:          0.0\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/lib/python2.7/dist-packages/ipykernel_launcher.py:6: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  \n",
      "/usr/local/lib/python2.7/dist-packages/ipykernel_launcher.py:36: RuntimeWarning: invalid value encountered in double_scalars\n"
     ]
    }
   ],
   "source": [
    "# test empty mask in ground truth\n",
    "mask_gt   = np.zeros((128,128,128), np.uint8)\n",
    "mask_pred = np.zeros((128,128,128), np.uint8)\n",
    "#mask_gt[50,60,70] = 1\n",
    "mask_pred[50,60,72] = 1\n",
    "surface_distances = compute_surface_distances(mask_gt, mask_pred, spacing_mm=(3,2,1))\n",
    "print(\"average surface distance: {} mm\".format(compute_average_surface_distance(surface_distances)))\n",
    "print(\"hausdorff (100%):         {} mm\".format(compute_robust_hausdorff(surface_distances, 100)))\n",
    "print(\"hausdorff (95%):          {} mm\".format(compute_robust_hausdorff(surface_distances, 95)))\n",
    "print(\"surface overlap at 1mm:   {}\".format(compute_surface_overlap_at_tolerance(surface_distances, 1)))\n",
    "print(\"surface dice at 1mm:      {}\".format(compute_surface_dice_at_tolerance(surface_distances, 1)))\n",
    "print(\"volumetric dice:          {}\".format(compute_dice_coefficient(mask_gt, mask_pred)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {
     "autoexec": {
      "startup": false,
      "wait_interval": 0
     },
     "base_uri": "https://localhost:8080/",
     "height": 238
    },
    "colab_type": "code",
    "executionInfo": {
     "elapsed": 476,
     "status": "ok",
     "timestamp": 1528999074739,
     "user": {
      "displayName": "Olaf Ronneberger",
      "photoUrl": "//lh5.googleusercontent.com/-Hp2EnPvUOmU/AAAAAAAAAAI/AAAAAAAAAFc/QR20zxrgwZ0/s50-c-k-no/photo.jpg",
      "userId": "110531354013919487725"
     },
     "user_tz": -60
    },
    "id": "bcgpOy5EqMLv",
    "outputId": "f929a69f-b1e9-4a51-f424-c02f7340aff7"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "average surface distance: (nan, nan) mm\n",
      "hausdorff (100%):         inf mm\n",
      "hausdorff (95%):          inf mm\n",
      "surface overlap at 1mm:   (nan, nan)\n",
      "surface dice at 1mm:      nan\n",
      "volumetric dice:          nan\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/lib/python2.7/dist-packages/ipykernel_launcher.py:6: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  \n",
      "/usr/local/lib/python2.7/dist-packages/ipykernel_launcher.py:7: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  import sys\n",
      "/usr/local/lib/python2.7/dist-packages/ipykernel_launcher.py:36: RuntimeWarning: invalid value encountered in double_scalars\n",
      "/usr/local/lib/python2.7/dist-packages/ipykernel_launcher.py:37: RuntimeWarning: invalid value encountered in double_scalars\n",
      "/usr/local/lib/python2.7/dist-packages/ipykernel_launcher.py:48: RuntimeWarning: invalid value encountered in double_scalars\n"
     ]
    }
   ],
   "source": [
    "# test both masks empty\n",
    "mask_gt   = np.zeros((128,128,128), np.uint8)\n",
    "mask_pred = np.zeros((128,128,128), np.uint8)\n",
    "#mask_gt[50,60,70] = 1\n",
    "#mask_pred[50,60,72] = 1\n",
    "surface_distances = compute_surface_distances(mask_gt, mask_pred, spacing_mm=(3,2,1))\n",
    "print(\"average surface distance: {} mm\".format(compute_average_surface_distance(surface_distances)))\n",
    "print(\"hausdorff (100%):         {} mm\".format(compute_robust_hausdorff(surface_distances, 100)))\n",
    "print(\"hausdorff (95%):          {} mm\".format(compute_robust_hausdorff(surface_distances, 95)))\n",
    "print(\"surface overlap at 1mm:   {}\".format(compute_surface_overlap_at_tolerance(surface_distances, 1)))\n",
    "print(\"surface dice at 1mm:      {}\".format(compute_surface_dice_at_tolerance(surface_distances, 1)))\n",
    "print(\"volumetric dice:          {}\".format(compute_dice_coefficient(mask_gt, mask_pred)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {
     "autoexec": {
      "startup": false,
      "wait_interval": 0
     }
    },
    "colab_type": "code",
    "id": "sVyRElyxYzzV"
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "default_view": {},
   "name": "Surface_distance_based_measures.ipynb",
   "provenance": [],
   "version": "0.3.2",
   "views": {}
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
